home *** CD-ROM | disk | FTP | other *** search
/ EnigmA Amiga Run 1999 March / EnigmA AMIGA RUN 35 (1999)(G.R. Edizioni)(IT)[!][issue 1999-03].iso / earcd / grafica / pvrgjpeg / leedct.c < prev    next >
C/C++ Source or Header  |  1999-01-01  |  9KB  |  469 lines
  1. /*************************************************************
  2. Copyright (C) 1990, 1991, 1993 Andy C. Hung, all rights reserved.
  3. PUBLIC DOMAIN LICENSE: Stanford University Portable Video Research
  4. Group. If you use this software, you agree to the following: This
  5. program package is purely experimental, and is licensed "as is".
  6. Permission is granted to use, modify, and distribute this program
  7. without charge for any purpose, provided this license/ disclaimer
  8. notice appears in the copies.  No warranty or maintenance is given,
  9. either expressed or implied.  In no event shall the author(s) be
  10. liable to you or a third party for any special, incidental,
  11. consequential, or other damages, arising out of the use or inability
  12. to use the program for any purpose (or the loss of data), even if we
  13. have been advised of such possibilities.  Any public reference or
  14. advertisement of this source code should refer to it as the Portable
  15. Video Research Group (PVRG) code, and not by any author(s) (or
  16. Stanford University) name.
  17. *************************************************************/
  18. /*
  19. ************************************************************
  20. leedct.c
  21.  
  22. This is the Byeong Gi Lee algorithm from IEEE Trans. Accoustics,
  23. Speech, and Signal Processing, Vol ASSP-32, No. 6, December 1984, pp.
  24. 1243 -1245.
  25.  
  26. ************************************************************
  27. */
  28.  
  29. /*LABEL leedct.c */
  30.  
  31. /*PUBLIC*/
  32.  
  33. extern void LeeIDct();
  34. extern void LeeDct();
  35.  
  36. /*PRIVATE*/
  37.  
  38. /* Standard Macros */
  39.  
  40. #define LS(r,s) ((r) << (s))
  41. #define RS(r,s) ((r) >> (s)) /* Caution with rounding... */
  42.  
  43. #define MSCALE(expr)  RS((expr),9)
  44.  
  45. #define IDCTSCALE(x) (((x<0) ? (x-8) : (x+8))/16);
  46. #define DCTSCALE(x) (((x<0) ? (x-8) : (x+8))/16);
  47.  
  48. /* Cos Table */
  49.  
  50. #define twoc1d4 724L
  51. #define twoc1d8 946L
  52. #define twoc3d8 392L
  53. #define twoc1d16 1004L
  54. #define twoc3d16 851L
  55. #define twoc5d16 569L
  56. #define twoc7d16 200L
  57. #define sqrt2 724L
  58.  
  59. /* 1/Cos Table */
  60.  
  61. #define itwoc1d4 362L
  62. #define itwoc1d8 277L
  63. #define itwoc3d8 669L
  64. #define itwoc1d16 261L
  65. #define itwoc3d16 308L
  66. #define itwoc5d16 461L
  67. #define itwoc7d16 1312L
  68. #define isqrt2 362L
  69.  
  70. #define x0 tx0
  71. #define x1 tx1
  72. #define x2 tx2
  73. #define x3 tx3
  74. #define x4 ex4
  75. #define x5 ex5
  76. #define x6 ex6
  77. #define x7 ex7
  78.  
  79. #define r0 rx0
  80. #define r1 rx1
  81. #define r2 rx2
  82. #define r3 rx3
  83.  
  84. #define s0 rx0
  85. #define s1 rx1
  86. #define s2 rx2
  87. #define s3 rx3
  88.  
  89. #define f0 ex0
  90. #define f1 ex1
  91. #define f2 ex2
  92. #define f3 ex3
  93.  
  94. #define g0 ex4
  95. #define g1 ex5
  96. #define g2 ex6
  97. #define g3 ex7
  98.  
  99. #define b1 gx0
  100. #define b2 gx0
  101. #define b3 gx1
  102.  
  103. #define a1 gx2
  104. #define a3 gx2
  105. #define c1 gx2
  106. #define c3 gx2
  107.  
  108. #define ihold gx1
  109.  
  110. /*START*/
  111.  
  112. /*BFUNC
  113.  
  114. LeeIDct is implemented according to the inverse dct flow diagram in
  115. the paper.  It takes two input arrays that must be defined before the
  116. call.
  117.  
  118. EFUNC*/
  119.  
  120. void LeeIDct(x,y)
  121.      int *x;
  122.      int *y;
  123. {
  124.   register int ex0,ex1,ex2,ex3,ex4,ex5,ex6,ex7;
  125.   register int tx0,tx1,tx2,tx3;
  126.   register int rx0,rx1,rx2,rx3;
  127.   register int gx0,gx1,gx2;
  128.   register int *iptr,*jptr;
  129.   register int i;
  130.  
  131.   /* Do rows */
  132.  
  133.   for(jptr=y,iptr=x,i=0;i<8;i++)
  134.     {
  135.       x0 = MSCALE(isqrt2*LS(*(iptr++),2));
  136.       x1 = LS(*(iptr++),2); 
  137.       x2 = LS(*(iptr++),2); 
  138.       x3 = LS(*(iptr++),2); 
  139.       x4 = LS(*(iptr++),2); 
  140.       x5 = LS(*(iptr++),2); 
  141.       x6 = LS(*(iptr++),2); 
  142.       x7 = LS(*(iptr++),2); 
  143.       
  144.       a1 = MSCALE(itwoc1d4*x4);
  145.       r0 = x0+a1;
  146.       r1 = x0-a1;
  147.       
  148.       a3 = MSCALE(itwoc1d4*(x2+x6));
  149.       r2 = MSCALE(itwoc1d8*(x2+a3));
  150.       r3 = MSCALE(itwoc3d8*(x2-a3));
  151.       
  152.       f0 = r0+r2;
  153.       f1 = r1+r3;
  154.       f2 = r0-r2;
  155.       f3 = r1-r3;
  156.       
  157.       b1 = x3+x5;
  158.       c1 = MSCALE(itwoc1d4*b1);
  159.       s0 = x1+c1;
  160.       s1 = x1-c1;
  161.       
  162.       b2 = x1+x3;
  163.       b3 = x5+x7;
  164.       c3 = MSCALE(itwoc1d4*(b2+b3));
  165.       s2 = MSCALE(itwoc1d8*(b2+c3));
  166.       s3 = MSCALE(itwoc3d8*(b2-c3));
  167.       
  168.       g0 = MSCALE(itwoc1d16*(s0+s2));
  169.       g1 = MSCALE(itwoc3d16*(s1+s3));
  170.       g2 = MSCALE(itwoc7d16*(s0-s2));
  171.       g3 = MSCALE(itwoc5d16*(s1-s3));
  172.       
  173.       *(jptr++) = f0+g0;
  174.       *(jptr++) = f1+g1;
  175.       *(jptr++) = f3+g3;
  176.       *(jptr++) = f2+g2;
  177.       
  178.       *(jptr++) = f2-g2;
  179.       *(jptr++) = f3-g3;
  180.       *(jptr++) = f1-g1;
  181.       *(jptr++) = f0-g0;
  182.     }
  183.  
  184.  
  185.   /* Do columns */
  186.  
  187.   for(i=0;i<8;i++)
  188.     {
  189.       jptr = iptr = y+i;
  190.  
  191.  
  192. #ifdef PVERSION
  193.  
  194.       x0 = MSCALE(isqrt2*(*(iptr)));
  195.       iptr += 8;
  196.       x1 = *(iptr); 
  197.       iptr += 8;
  198.       x2 = *(iptr); 
  199.       iptr += 8;
  200.       x3 = *(iptr); 
  201.       iptr += 8;
  202.       x4 = *(iptr); 
  203.       iptr += 8;
  204.       x5 = *(iptr); 
  205.       iptr += 8;
  206.       x6 = *(iptr); 
  207.       iptr += 8;
  208.       x7 = *(iptr); 
  209.  
  210. #else
  211.  
  212. #undef x1
  213. #undef x2
  214. #undef x3
  215. #undef x4
  216. #undef x5
  217. #undef x6
  218. #undef x7
  219.  
  220. #define x1 iptr[8]
  221. #define x2 iptr[16]
  222. #define x3 iptr[24]
  223. #define x4 iptr[32]
  224. #define x5 iptr[40]
  225. #define x6 iptr[48]
  226. #define x7 iptr[56]
  227.  
  228.       x0 = MSCALE(isqrt2*(*iptr));
  229.  
  230. #endif
  231.       
  232.       a1 = MSCALE(itwoc1d4*x4);
  233.       r0 = x0+a1;
  234.       r1 = x0-a1;
  235.       
  236.       a3 = MSCALE(itwoc1d4*(x2+x6));
  237.       r2 = MSCALE(itwoc1d8*(x2+a3));
  238.       r3 = MSCALE(itwoc3d8*(x2-a3));
  239.       
  240.       f0 = r0+r2;
  241.       f1 = r1+r3;
  242.       f2 = r0-r2;
  243.       f3 = r1-r3;
  244.       
  245.       b1 = x3+x5;
  246.       c1 = MSCALE(itwoc1d4*b1);
  247.       s0 = x1+c1;
  248.       s1 = x1-c1;
  249.       
  250.       b2 = x1+x3;
  251.       b3 = x5+x7;
  252.       c3 = MSCALE(itwoc1d4*(b2+b3));
  253.       s2 = MSCALE(itwoc1d8*(b2+c3));
  254.       s3 = MSCALE(itwoc3d8*(b2-c3));
  255.       
  256.       g0 = MSCALE(itwoc1d16*(s0+s2));
  257.       g1 = MSCALE(itwoc3d16*(s1+s3));
  258.       g2 = MSCALE(itwoc7d16*(s0-s2));
  259.       g3 = MSCALE(itwoc5d16*(s1-s3));
  260.  
  261.       ihold = f0+g0;
  262.       (*jptr) = IDCTSCALE(ihold);
  263.       jptr += 8;
  264.       ihold = f1+g1;
  265.       (*jptr) = IDCTSCALE(ihold);
  266.       jptr += 8;
  267.       ihold = f3+g3;
  268.       (*jptr) = IDCTSCALE(ihold);
  269.       jptr += 8;
  270.       ihold = f2+g2;
  271.       (*jptr) = IDCTSCALE(ihold);
  272.       jptr += 8;
  273.       ihold = f2-g2;
  274.       (*jptr) = IDCTSCALE(ihold);
  275.       jptr += 8;
  276.       ihold = f3-g3;
  277.       (*jptr) = IDCTSCALE(ihold);
  278.       jptr += 8;
  279.       ihold = f1-g1;
  280.       (*jptr) = IDCTSCALE(ihold);
  281.       jptr += 8;
  282.       ihold = f0-g0;
  283.       (*jptr) = IDCTSCALE(ihold);
  284.  
  285.     }
  286. }
  287.  
  288. #undef f0
  289. #undef f1
  290. #undef f2
  291. #undef f3
  292.  
  293. #undef g0
  294. #undef g1
  295. #undef g2
  296. #undef g3
  297.  
  298. #undef r0
  299. #undef r1
  300. #undef r2
  301. #undef r3
  302.  
  303. #undef s0
  304. #undef s1
  305. #undef s2
  306. #undef s3
  307.  
  308. #define f0 rx0
  309. #define f1 rx1
  310. #define f2 rx2
  311. #define f3 rx3
  312.  
  313. #define r0 rx4
  314. #define r1 rx5
  315. #define r2 rx6
  316. #define r3 rx7
  317.  
  318. #define g0 sx0
  319. #define g1 sx1
  320. #define g2 sx2
  321. #define g3 sx3
  322.  
  323. #define s0 rx4
  324. #define s1 rx5
  325. #define s2 rx6
  326. #define s3 rx7
  327.   
  328.  
  329. /*BFUNC
  330.  
  331. LeeDct is implemented by reversing the arrows in the inverse dct flow
  332. diagram.  It takes two input arrays that must be defined before the
  333. call.
  334.  
  335. EFUNC*/
  336.  
  337. void LeeDct(x,y)
  338.      int *x;
  339.      int *y;
  340. {
  341.   register int rx0,rx1,rx2,rx3,rx4,rx5,rx6,rx7;
  342.   register int sx0,sx1,sx2,sx3;
  343.   register int hold,c2;
  344.   register int *iptr,*jptr;
  345.  
  346. #undef x0
  347. #undef x1
  348. #undef x2
  349. #undef x3
  350. #undef x4
  351. #undef x5
  352. #undef x6
  353. #undef x7
  354.  
  355. #define x0 iptr[0]
  356. #define x1 iptr[1]
  357. #define x2 iptr[2]
  358. #define x3 iptr[3]
  359. #define x4 iptr[4]
  360. #define x5 iptr[5]
  361. #define x6 iptr[6]
  362. #define x7 iptr[7]
  363.  
  364.   for(jptr=y,iptr=x;iptr<x+64;jptr+=8,iptr+=8)
  365.     {
  366.       f0 = x0+x7;
  367.       g0 = MSCALE(twoc1d16*(x0-x7));
  368.       f1 = x1+x6;
  369.       g1 = MSCALE(twoc3d16*(x1-x6));
  370.       f3 = x2+x5;
  371.       g3 = MSCALE(twoc5d16*(x2-x5));
  372.       f2 = x3+x4;
  373.       g2 = MSCALE(twoc7d16*(x3-x4));
  374.  
  375.       r0 = f0+f2;
  376.       r1 = f1+f3;
  377.       r2 = MSCALE(twoc1d8*(f0-f2));
  378.       r3 = MSCALE(twoc3d8*(f1-f3));
  379.  
  380.       jptr[0] = MSCALE(sqrt2*(r0+r1));
  381.       jptr[4] = MSCALE(twoc1d4*(r0-r1));
  382.       jptr[2] = hold = r2+r3;
  383.       jptr[6] = MSCALE(twoc1d4*(r2-r3))-hold;
  384.  
  385.       s0 = g0+g2;
  386.       s1 = g1+g3;
  387.       s2 = MSCALE(twoc1d8*(g0-g2));
  388.       s3 = MSCALE(twoc3d8*(g1-g3));
  389.       
  390.       jptr[1] = hold = s0+s1;
  391.       c2 = s2+s3;
  392.       jptr[3] = hold = c2-hold;
  393.       jptr[5] = hold = MSCALE(twoc1d4*(s0-s1))-hold;
  394.       jptr[7] = MSCALE(twoc1d4*(s2-s3))-c2-hold;
  395.     }
  396.  
  397. #undef x0
  398. #undef x1
  399. #undef x2
  400. #undef x3
  401. #undef x4
  402. #undef x5
  403. #undef x6
  404. #undef x7
  405.  
  406. #define x0 iptr[0]
  407. #define x1 iptr[8]
  408. #define x2 iptr[16]
  409. #define x3 iptr[24]
  410. #define x4 iptr[32]
  411. #define x5 iptr[40]
  412. #define x6 iptr[48]
  413. #define x7 iptr[56]
  414.  
  415. #define y0 rx0
  416. #define y1 rx0
  417. #define y2 rx1
  418. #define y3 rx1
  419. #define y4 rx2
  420. #define y5 rx2
  421. #define y6 rx3
  422. #define y7 rx3
  423.  
  424.   for(jptr=y,iptr=y;iptr<y+8;jptr++,iptr++)
  425.     {
  426.       f0 = x0+x7;
  427.       g0 = MSCALE(twoc1d16*(x0-x7));
  428.       f1 = x1+x6;
  429.       g1 = MSCALE(twoc3d16*(x1-x6));
  430.       f3 = x2+x5;
  431.       g3 = MSCALE(twoc5d16*(x2-x5));
  432.       f2 = x3+x4;
  433.       g2 = MSCALE(twoc7d16*(x3-x4));
  434.  
  435.       r0 = f0+f2;
  436.       r1 = f1+f3;
  437.       r2 = MSCALE(twoc1d8*(f0-f2));
  438.       r3 = MSCALE(twoc3d8*(f1-f3));
  439.       
  440.       y0 = MSCALE(sqrt2*(r0+r1));
  441.       y4 = MSCALE(twoc1d4*(r0-r1));
  442.       y2 = r2+r3;
  443.       y6 = MSCALE(twoc1d4*(r2-r3))-y2;
  444.  
  445.       jptr[0] = DCTSCALE(y0);
  446.       jptr[16] = DCTSCALE(y2);
  447.       jptr[32] = DCTSCALE(y4);
  448.       jptr[48] = DCTSCALE(y6);
  449.       
  450.       s0 = g0+g2;
  451.       s1 = g1+g3;
  452.       s2 = MSCALE(twoc1d8*(g0-g2));
  453.       s3 = MSCALE(twoc3d8*(g1-g3));
  454.       
  455.       y1 = s0+s1;
  456.       c2 = s2+s3;
  457.       y3 = c2-y1;
  458.       y5 = MSCALE(twoc1d4*(s0-s1))-y3;
  459.       y7 = MSCALE(twoc1d4*(s2-s3))-c2-y5;
  460.  
  461.       jptr[8] = DCTSCALE(y1);
  462.       jptr[24] = DCTSCALE(y3);
  463.       jptr[40] = DCTSCALE(y5);
  464.       jptr[56] = DCTSCALE(y7);
  465.     }
  466. }
  467.  
  468. /*END*/
  469.